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Novel electronic order is found theoretically for a system where even number of localized 
electrons per site are coupled with conduction electrons. For precise quantitative study, a variant 
of the Kondo lattice model is taken with crystalline electric field (CEF) singlet and triplet 
states for each site. Using the dynamical mean-field theory combined with the continuous-time 
quantum Monte Carlo method, a staggered order with alternating Kondo and CEF singlets is 
identified for a case with one conduction electron per site being distributed in two conduction 
bands each of which is quarter-filled. This electronic order accompanies a charge density wave 
(CDW) of conduction electrons that accumulate more on Kondo-singlet sites than on CEF- 
singlet sites. Possible relevance of the present order to the scalar order in PrFe4Pi2 is discussed. 
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1. Introduction 

Dichotomy between itinerant and localized characters 
of /-electrons leads to intriguing phenomena in heavy 
fermion systems. With f 1 configuration as in Ce com- 
pounds, the behavior of the system is determined by the 
competition between the Kondo effect and the RKKY in- 
teraction. The former gives itinerancy to /-electrons. If / 
electrons are localized, the RKKY interaction gives rise 
to a magnetic order. Depending on the relative strength 
of the Kondo effect and the RKKY interaction, the 
ground state can either be magnetic or nonmagnetic. 
This picture was first spelled out by Doniach, 1 ' and is 
now recognized well. On the other hand, in systems with 
two /-electrons per site as in some Pr and U compounds, 
the distinction between itinerant and localized characters 
is much more subtle. It remains to be clarified how the 
dichotomy plays a role in rich behaviors found in real 
systems such as PrFe4Pi2 and URu2Si2- 2 ' 3 ' 

A typical situation in / 2 systems is the case where 
the ground state is the crystalline electric field (CEF) 
singlet. If /-electrons interact strongly with conduction 
electrons, the ground state can be a collective Kondo 
singlet involving higher CEF levels. Also in this case, 
/-electrons acquire itinerancy because of the Kondo ef- 
fect. Competition between these two singlets, namely di- 
chotomy between itinerant and localized characters of 
/-electrons, may give rise to rich physics with exotic or- 
dered phases. Theoretically, there are many discussions 
in impurity systems with the CEF singlet. 4 ~ 9 ' In order 
to discuss electronic orders, however, investigation of pe- 
riodic systems is necessary. 

Let us consider the case where CEF states form a 
quasi-quartet composed of a singlet ground state and a 
first-excited triplet. We represent the singlet-triplet lev- 
els at each site i in terms of two pseudo spins S-yi with 
7 = 1,2. Correspondingly, we introduce two conduction 
bands 7=1,2 where electrons interact with pseudo spins 
with the same orbital index 7 by the exchange interaction 

*E-mail: hoshino@cmpt.phys. tohoku. ac.jp 



J 7 > 0. Then the model reads 

hycr 27 



07* 



+ A^2S U -S 2 



(1) 



where s cl i denotes the spin of conduction electrons at 
site i. The CEF splitting A is simulated by the coupling 
between Su and S 2 i as shown in the third term . We 
call this model the "two-band singlet-triplet Kondo lat- 
tice model" (2BSTKLM) in this paper. This model at 
half filling has been investigated in one dimension by the 
density- matrix renormalization group. 10 -' 

In the case of A = 0, the 2BSTKLM corresponds to the 
two independent Kondo lattice models. 11,12 ' Recently, a 
charge density wave (CDW) transition has been found in 
the Kondo lattice model at quarter filling. 13 ' However, 
this ordered state cannot be the ground state since local- 
ized spins keep finite entropy even below the transition 
temperature. In this paper, we investigate an effect of the 
CEF splitting on the CDW state. We shall demonstrate 
that the CEF splitting stabilizes a non-magnetic order- 
ing, and leads to a novel ground state, which we call the 
"staggered Kondo-CEF singlet order" in the following. 

In order to investigate the ordering, we employ the 
dynamical mean- field theory (DMFT), 14 ' which can be 
extended for the ordered state in the Kondo lattice sys- 
tems. 15 ' 16 ' The DMFT takes full account of on-site cor- 
relations, and becomes exact in infinite dimensions. Since 
both Kondo and CEF effects are dominantly local, the 
DMFT is suitable for our purpose. As the impurity solver 
associated with the DMFT, we use the continuous-time 
quantum Monte Carlo method (CT-QMC). 17 " 19 ' 

In §2 we describe the formalism how to derive suscep- 
tibilities in the ordered phase with A and B sublattices 
in the framework of the DMFT . §3 is devoted to techni- 
cal aspects about the CT-QMC. The main results of this 
paper is given in §4 and §5, where we derive the phase 
diagram, and some physical quantities such as suscep- 
tibilities, order parameters, local correlation functions, 
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and renormalized density of states. We discuss in §6 how 
the characteristics of the system is compared with the 
ordinary Kondo lattice, and how the present results are 
relevant to understanding the scalar order in PrFe4Pi2. 
The summary of the paper is given in §7, and some tech- 
nical details related to §2 are given in Appendix. 

2. Susceptibilities with Two-Sublattices 

An instability toward ordered phase is signalled by di- 
vergence of the corresponding susceptibility. In this sec- 
tion, we generalized the DMFT formalism 14 ) so that we 
can investigate the instability including two-sublatticc 
systems. We first consider the simplest case of non- 
interacting conduction electrons with nearest-neighbor 
hopping in a bipartite lattice, where the condition 
£k+Q = —£k with Q — (tt, 7r, • • • ) is satisfied. We use 
the suffix a for spin and orbital indices. In particular, 
a = (7,0-) in the 2BSTKLM given by eq. (1). We in- 
troduce sublattice annihilation operators CkA(c k B) for 
A(B)-sublattice with wave vector k. The first term in 
cq. (1) is rewritten in terms of the sublattice operators 
as 



Ho 



k ex 



CfcBc 



A 

-fcBa 



CfeAc 



CfcAa + c fcBa C kB 



(2) 



where the summation ^ is taken over wave vectors k 
belonging to the reduced Brillouin zone with sublattices, 
i.e., half of the original Brillouin zone. 

The effect of terms other than Hq in eq.(l) appears 
as the self-energy £> Q (z) of conduction electrons, where 
A = A, B indicates the label for the sublattice. Note that 
the self energy is spatially local in the DMFT. In terms 
of the quantity Cxa(z) = z + [i — T,\ a (z), the full Green 
function is given by 

Gfca ( z ) = 

- [^(z)6 xx ,+e ha (l-6xx>)]. (3) 



where A indicates the complementary component such as 
A = B. 

Let us assume a symmetric density of states: p(e) = 
(N/2)- 1 J2'k S ( £ - £ k) = P(s). Note that p(e) is the 
same as the density of states in the original Brillouin 
zone. In the symmetric case, the local Green function 
becomes diagonal with respect to the sublattice label. 
This is easily seen from eq. (3) where the off-diagonal part 
is an odd function of s kai and vanishes by summation 
over k. Then the local Green function is given by 



(4) 



Magnetic and charge susceptibilities can be evaluated 
from the two-particle Green function, which is given in 
the imaginary time domain by 



Xaa'( r l,T 2 ,T3,T 4 ) 



N/2 ^ 
1 fc.fc 



\_{ T rCl Xa {Tl)c kXa (T 2 )4, A , Q , (T 3 ) Ck , X , a , (T4)) 




Xiai 



Vet' 



Fig. 1. Bethe-Salpeter equation for the two-particle Green func- 
tion in the two-sublattice formulation. Summation should be 
taken over the sublattice index Ai and the spin-orbital index 



-( T rcj c Aa( T l) C fcAQ(T 2 ))(rr4'A' Q '( T 3) C fc'A'a'(T4)) ■ (5) 

Here we only consider the uniform component in the re- 
duced Brillouin zone. The uniform and staggered com- 
ponents in the original Brillouin zone can be calculated 
by the linear combinations as 



,unif ,sta 



X 



1(X AA + X BB ±2 X AB ). (6) 



The Fourier transform for the two-particle Green func- 
tion (5) is defined by 

1 ff> 

X(i£ n ,i£ n /;ii/ m ) = -J / dTidT 2 dT 3 dT4 
P Jo 



xe 



i£n(T2— ti) ie„/(r 4 -T 3 ) (r 2 -r 3 ) 



x(ti,t 2 ,T3,t 4 ), (7) 



where e n and v m are the fermionic and bosonic Mat- 
subara frequencies, respectively. The susceptibility which 
represents a response to external fields is derived from 
X(ie ni ie n >;iv m ) as 



,' fern ) = 4 *W fen fen'', 'Wm ) ■ 



(8) 



We use the Bethe-Salpeter equation to relate Xua> to 
the effective impurity model. Noting that the vertex part 
is local and now depends on the sublattice, we obtain the 
following equation: 

Xaa'(i£n,i£n',hVm) = Xa^ fen, '^m)S aa ' $nn> 

+ E X Q ot XXl fenAv m ) 
\\ct\ni 

XTaax fen, '^m ! '^m)Xa\a' ' ie »' 5 {v ™)- (9) 

The function x° is the two-particle Green function with- 
out the vertex part, and defined by 

X ° a XX ' fe n ;iv m ) = -J^J2' G kafen)G X ka X fe n +^ m ). 

(10) 

Equation (9) is graphically shown in Fig. 1. 

We evaluate the local vertex T A using the effective 
impurity model. We define the local two-particle Green 
function as 

Xloc,aQ'( r l, T 2 ,T 3 ,T 4 ) = 
( r rcL( T l) C Aa(T2)4 Q ,(T 3 )cA Q '(T4)) 

- ( T rci Q (Ti)cA Q (T 2 ))(T T 4 Q ,(r 3 )cAa'(T4)), (11) 

where c\ a = (7V/2) -1 / 2 J2' k c k\a is the annihilation oper- 
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ator at the origin site for each sublattice. In the effective 
impurity model, the Bethe-Salpeter equation for Xio C aa' 
is expressed by local quantities as 

+ zZ X°oc,a( i£ «>™) 
aini 

X r aa 1 ( i£ n. i£ n 1 ; it 'm)Xloc,a 1 a'( i£ ni) i £n';i^m), (12) 

where 

X°oc,a( i£ ™; {Vm ^> = - G ^oc,a(i £ ™)Gioc,a(i £ n + iv m ). (13) 

We note that, in general, the local two-particle Green 
function in cq. (11) has the off-diagonal element such 
as Xkic- However with the symmetric property p(e) = 
p(—e), the off-diagonal elements vanish because is 
diagonal. 

In the DMFT, the vertices in eqs. (9) and (12) are 
the same. Hence, we can obtain x AA by solving these 
equations simultaneously. To make the notation simple, 
we use a matrix form with respect to (n, a). Then eq. 
(12) is rewritten as 



0,Ai-l 
loci 



IxLr = \x 

Similarly, eq. (9) is rewritten as 



(14) 



7 AA 
BA 



X 



X ABX-1 /0,AA X 0,AB 
^BB I — I X°' BA X°' BB 



(15) 

Consequently, the problem is reduced to the calculation 
of the local two-particle Green function Xk, c an d Xk> c 
in the effective impurity model. These can be evaluated 
by the two-particle i-matrix as discussed in ref. 20. The 
details of the calculations of x° an d X? oc are g iven m 
Appendix. 

3. CT-QMC Algorithm 

In the DMFT, a periodic model is mapped to an ef- 
fective impurity model. 14 ' Here we explain how to apply 
the CT-QMC to the present impurity model. An impu- 
rity version of the 2BSTKLM (1) is written as follows: 

(16) 
(17) 

(18) 



T-L — 7i c + Hf + "Hint, 

He = ^ £ ky C kja C k'Y<T + ^ 



ky<7 



7(7 



H f = -AP s + -, 



■yaa' 7 

(19) 

where v 7 = — J 7 /4, and Aj CT , = \"fa){"fa'\ is the X- 
operator for the localized states, and P s = — Si ■ S2 + 1/4 
is the projection operator onto the CEF singlet state. The 
parameter a 7 is chosen as for ferromagnetic interaction 
(J 7 < 0) and 1 for antiferromagnetic coupling (J 7 > 0) 
in order to avoid the minus sign problem as noted in ref. 
9. The constant term may be neglected in the simulation. 
The partition function Z is expanded with respect to 



"Hint using the formula 
Z = Tr \ T T e~ m ° exp 



- j\rH\ nt {r) j, (20) 



where the non-interacting Hamiltonian Hq is defined 
by %o = 7i c + Hf. The suffix T denotes the interac- 
tion picture: W? nt (r) = e TUa n- m t^ rUa ■ The CT-QMC 
evaluates the perturbation expansion by Monte Carlo 
method. 17 ' 18 ^ The trace over conduction electrons is 
computed with the aid of Wick's theorem. For the local- 
ized spins, on the other hand, we evaluate the following 
expression: 



W f = (Xjf a ,(r k )...X^,(n)) f 



71 1 



(21) 



Here k is the perturbation order, and we consider the 
imaginary-time set with f3 > T}. > ■ ■ ■ > T\ > 0. The 
average is defined by (• • • ) f — Trj[e _/3W/ • • - ]/Zf with 
Zf = 3 + e ,9A . The efficient algorithm for the Kondo 
model 9 ' 19 ) using the "segment" cannot be applied to the 
present model, because the localized state \^a) is not 
an eigenstate of Hf. Therefore, we must multiply the k 
matrices to evaluate eq. (21). 
We rewrite eq. (21) as 



— Tr f [F k F k ^i ■ ■ ■ Fi] , 



W f 

Fi = P {n +1 n)X^, 



(22) 
(23) 



where we have introduced p(r) = e~ rH f , and r k +i is de- 
fined by T k +i = (3 + Ti. This multiplication process takes 
the main part of the computational time, and the cal- 
culation becomes heavy compared to the method using 
segments. This difficulty can be somewhat reduced by us- 
ing the so-called tree algorithm or binning algorithm. 21 ) 
Although the direct multiplication takes the computing 
time of 0{k), the tree algorithm and binning algorithm 
give O(logfc) and 0(y/k), respectively. We have imple- 
mented the binning algorithm, since it is much simpler 
than the other. 

For evaluating eq. (22), it is convenient to choose the 
basis that diagonalizes p{r). We denote this eigenstates 
by (|s), |t+), I tO), |t— )), and p{r) becomes 



P(r) = 



V 



V 



(24) 



U = 



(25) 



Correspondingly, the A-operator Aj CT , expressed in the 
basis (| tt), I U), I It), I U» is replaced by UXj a ,U-\ 
where 

fo 1/V2 -1/V2 0\ 

10 

1/V2 l/y/2 

Vo 01/ 

Using this representation, it can be shown that negative 
weight does not appear up to the second order of J 7 . 
The absence in fact persists to higher orders empirically. 
In the presence of magnetic fields, however, the negative 
sign may arise. 
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Fig. 2. Equal-time correlation (Si ■ s c i) as a function of temper- 
ature. The parameters are Ji = J2 = A = 0.5, 1, 2, 4, 8, 12 
from top to bottom. The dotted line shows the value (Si ■ s c i) — 
—0.683 for the ground state in the strong coupling limit . 



Next we discuss how to calculate physical quantities. 
The algorithm for conduction electrons is the same as 
the previous CT-QMC methods. For localized part, the 
quantity (T t A(t)B) with A and B being operators for 
localized states is given by the following expression: 

(TrA(r)B) = 

±- (Tr f [F k ■ ■ ■ F i+1 A\r r i+1 )F ; • • • F.B] > mc , (26) 

where (• • • )mc means the Monte Carlo average, and we 
have assumed r £ [Tj,Tj + i]. With use of this formula, 
we can evaluate quantities such as (Si ■ S2} and time- 
dependent correlation functions. We note that there is 
another method to calculate the correlation functions us- 
ing the inverse matrix of the Green function for conduc- 
tion electrons. 19 -* 

The correlation between localized and conduction 
spins can be evaluated by using another method. We 
write the Hamiltonian as % = Ho + A"Hi nt , and the 
partition function as Z\. With the expansion Z\ = 
J2kLo ^ k Zk, we obtain the following equation: 

A A=l k=0 

From eq. (20), on the other hand, the left-hand side of 

(27) is given by —(3 (Hint). Thus we derive the formula 

(U int ) = -(k) MC //3. (28) 

Since the Hamiltonian Hint includes the interaction be- 
tween conduction and localized spins, we can evaluate 
(Sj ■ s C7 ) from this expression. Note that the expression 

(28) cannot be applied to the time-dependent correlation 
functions. 

At the end of this section, we show an exemplary re- 
sult of our simulation to check the accuracy. Figure 2 
shows the equal-time correlation (Si • s c i) as a function 
of temperature with Ji = J 2 = A. Here we have used the 
rectangular density of states po(e) = 6(D — \e\)/2D, and 
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put D = 1. In the strong-coupling limit J 7 , A ^ D, the 
impurity model given by eq. (16) is reduced to four-spin 
system since the kinetic energy term can be neglected. 
Diagonalizing this Hamiltonian, we obtain (Si • s c \) — 
-(1 + \/3)/4 ~ -0.683 at T = 0. The result in Fig. 2 
tends to this value, and (Si ■ S2) tends to —1/4. Even 
though the CT-QMC is based on the expansion from the 
weak coupling, we have confirmed that it can reproduce 
the strong coupling limit. 

4. Phase Diagram 

In the rest of this paper, we put Ski = £k2 = £k, and 
use a tight-binding band on a hypercubic lattice. The 
density of states is given by 



Po(e) = 




with D — 1 as a unit of energy. This band has a perfect 
nesting property with the wave vector Q at half filling. 14 ) 
In order to discuss the effect of the CEF splitting, we fix 
the strength of the interaction as J\ = J 2 = J = 0.8 
and vary the CEF splitting A. Note that the definition 
of J is different from that in rcf. 13 by factor 2. We 
also fix the number of conduction electrons per site as 
n c = 1, which corresponds to the quarter filling of both 
conduction bands. 

4-1 Divergence of Susceptibilities 

To determine the phase diagram, we search for insta- 
bility of normal states in terms of divergent response of 
conduction electrons. The magnetic and charge response 
functions are defined by 

Xc, s = ^(xn + xw ±2 xn)> ( 30 ) 

where the indices 'c' and 's' indicate charge and spin 
channel, respectively. In a similar manner, we define '+' 
and ' — ' channels from the orbital index 

X± = 7j(Xll+X22±2xi 2 )- ( 31 ) 

Combining with the spatial dependence defined in 
eq. (6), we consider the following eight susceptibilities: 
x3: n c if ,X± n s if ,xS S ,X±s aS - Figure 3 illustrates the elec- 
tronic orders in the strong coupling limit, each of which 
is probed by the corresponding susceptibility. 

Figure 4 shows the temperature dependence of inverse 
susceptibilities for A = and A = 0.1 . Note that the 
case with A = shown in Fig. 4(a) corresponds to the 
pair of ordinary Kondo lattice models at quarter fill- 
ing. In this case, \+ an d X- are equivalent because of 
X12 = 0. For comparison, we include in Fig. 4(a) the high 
temperature form of the susceptibility given by 

Xhigh^W— ^-77-) = /3/(Mo)/(-Mo), (32) 

where f(x) = l/(e l3x + 1) is the Fermi distribution func- 
tion and po is the chemical potential without interaction. 
The susceptibilities of the 2BSTKLM show good agree- 
ment with Xhigh in the high temperature region. 

In Fig. 4(a), x±c S tends to diverge at T = T CDW ~ 
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(a) unif, +c 
1 



(e) stag, +c 
1 



(b) unif, -c 
1 



(f) stag, -c 
1- 



(c) unif, +s 



(g) stag, +s 



rrm:nm 

(d) unif, -s (h) stag, -s 

imrr :trm 



Fig. 3. Illustration of electronic orders probed by divergence of 
susceptibilities. The labels 1 and 2 are the orbital indices of con- 
duction electrons. 



0.015, which indicates the CDW instability. In this or- 
dered phase, two states illustrated in Figs. 3(e) and 
(f) are degenerate. This degeneracy is resolved by the 
singlet-triplet splitting. The CEF splitting suppresses 
X S -^ S an( i enhances x+c S leading to a higher transition 
temperature. Figure 4(b) demonstrates this situation, 
giving Xcdw ~ 0.028. As will be shown in §5.1, this 
CDW state is identified as the staggered Kondo-CEF 
singlet order. Although x™ lf also diverges at T ~ 0.012, 
this is not meaningful because the ordered state probed 
by x™ lf should be the reference state below Tcdw- 

Let us now focus on x+c S tna t nas been found to di- 
verge. Figure 5(a) shows the temperature dependence of 
X+c g f° r various values of A. The transition temperature 
Tcdw increases with increasing A. For A > 0.3, on the 
other hand, it is more convenient to consider the sus- 
ceptibility as a function of A with fixed T because of the 
sharp decrease of Tcdw- Fig. 5(b) shows that the criti- 
cal value of A increase with decreasing temperature. As 
clearly seen from this figure, the transition points have 
almost no difference between T = 0.015 and T = 0.01. 
Then we conclude that there is no ordering for A > 0.41. 

In this way, we obtain the T-A phase diagram as shown 
in Fig. 6 for J — 0.8 with one conduction electron per 
site. The increase of Tcdw in the small A region indi- 
cates a stabilization of the staggered Kondo-CEF singlet 
order by the CEF splitting. Sufficiently large A destroys 
the order, and the system becomes normal state with 
the CEF singlet. The transition changes from the second 
order to the first order around A = 0.4, which will be dis- 
cussed in the next subsection. For sufficiently small A, 
we expect a magnetically ordered ground state, as in the 
ordinary Kondo lattice, because of the residual entropy 
of the localized spins. This region is indicated qualita- 
tively in Fig. 6 as "Magnetic". More details about this 
aspect will be discussed in §6.1. We have also taken other 
values of coupling constant such as J = 0.6 or 1.0, and 
found no qualitative change in the overall behavior. 




unif,±,c 
unif,±,s 
stag,±,c 
stag,±,s 
high temp, limit 




Fig. 4. (color online) Inverse susceptibilities for (a) A = and (b) 
A = 0.1. In (a), + and — channels are equivalent since there is 
no off-diagonal component Xi2- The high-temperature behavior 
is given in eq. (32). 



4-2 Growth of Order Parameter and Hysteresis 

We discuss the order parameter defined by ti c q = 
(n c A — n C B)/2. Figure 7 shows the temperature de- 
pendence of n c Q. The transition temperature TcdWi 
which is determined by divergence of the susceptibility, 
is also shown in Fig. 7. The order parameter grows as 
(Tcdw ~ T) 1 / 2 as in the usual mean-field theory. We 
note that the initial rise around the transition tempera- 
ture becomes sharper with increasing A. 

Next we show the presence of first-order transition by 
deriving the order parameter as a function of A at fixed 
temperature. Figure 8 shows the result at T = 0.01. In 
the region with small A, n C Q becomes larger for larger 
CEF splitting, which indicates stabilization of the stag- 
gered Kondo-CEF singlet order. This behavior corre- 
sponds to the increase of Tcdw shown in Fig. 6. The 
magnitude of ti c q becomes the maximum at A ~ 0.25, 
and slowly decreases with larger splitting. 

As shown in Fig. 8, we have observed the hysteresis 
around A ~ 0.42. This is a characteristic for the first- 
order transition. When approaching from smaller A, we 
have used the initial condition with a staggered chemical 
potential for the effective medium of the DMFT. On the 
other hand, when approaching from larger CEF splitting, 
common chemical potential is used for both sublattices. 
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0.04 



J = 0.8 



Fig. 5. (color online) Inverse susceptibility as functions 

of (a) temperature T and (b) CEF splitting A. 



The shaded area in Fig. 6 shows the region of hysteresis. 
Note that the hysteresis is observed only in the region 
with T < 0.022. This means the existence of the critical 
point from second- to first-order transitions. More details 
are discussed in terms of the susceptibility in the next 
subsection. 

4-3 Critical Point from First- to Second- Order Transi- 
tions 

We study the hysteresis in more detail near A ~ 0.4. 
Figure 9 shows the inverse susceptibility at T = 0.015 in 
both disordered and ordered phases. For calculation of 
the susceptibility in the ordered phase, the formulation 
given in §2 is used. If the transition is of second order, 
both susceptibilities must diverge at the same transition 
point . The result in Fig. 9 shows divergent susceptibility 
in the normal state at A = Acdw, while the staggered 
Kondo-CEF singlet order phase has the finite susceptibil- 
ity at this point. Hence this difference clearly shows the 
first-order nature of the transition. It is difficult to reach 
the divergence of the susceptibility from the staggered 
Kondo-CEF singlet order phase since tiny statistical er- 
rors destroy the meta-stable ordered states. 

Now let us derive the critical point where the transition 
changes from second to first order. Following the proce- 
dure similar to that shown in Fig. 9, we derive Acdw 
for a given temperature by extrapolating 1 /x+c g to zero 




Fig. 6. Phase diagram of the 2BSTKLM with J = 0.8 and one 
conduction electron per site. The phase labelled "CDW" is iden- 
tified as the staggered Kondo-CEF singlet order. The shaded area 
shows the region of first-order transition, which ends at the crit- 
ical point. See text for details of derivation. Magnetic ordering 
is expected with small A and T as indicated by "Magnetic" . 
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Fig. 7. (color online) Temperature dependence of the order pa- 
rameter n c Q (left scale) as a function of temperature. Also shown 
is the inverse susceptibility [x+c S ] _1 (right scale). 



Of 




Fig. 8. Order parameter vs A with fixed temperature at T = 
0.01. The arrows show the hysteresis that occurs depending on 
whether we approach from smaller or larger A. 
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0.36 0.38 0.4 0.42 0.44 0.46 0.48 0.5 
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Fig. 9. Inverse susceptibility with temperature fixed at T = 0.015 
both inside and outside the ordered phase. With decreasing A 
in the disordered phase, l/x+c S tends to zero at A = Acdw> 
which is obtained by extrapolation. While 1 /x+c S derived in the 
ordered phase is finite at A = Aqdw ■ 
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Fig. 11. (color online) Local magnetic susceptibilities Xm vs 
temperature. With A = 0, the system decomposes into two in- 
dependent Kondo lattices, and the off-diagonal component Xm 
vanishes. The labels A and B show the electron-rich and electron- 
poor sites, respectively. 
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Fig. 10. Inverse susceptibility in the ordered phase as a function 
of temperature. The value A at each temperature is taken to be 
Aqdw as explained in the text. The dotted line is guided to the 
eye. 



from large-A side. Then we calculate x+c S at A = Acdw 
in the ordered phase taking the same temperature. Fig- 
ure 10 shows the result of x+c S calculated in this way 
for different values of temperature. The apparent scat- 
tering of the data for T < 0.018 comes mainly from 
ambiguity of extrapolation of the inverse susceptibility. 
With increasing temperature, the susceptibility becomes 
larger, and tends to diverge at T = T cr ~ 0.021 with 
A = A cr ~ 0.39. Since the susceptibility in the ordered 
phase must diverge also at Acdw if the transition is of 
second-order, we conclude that T cr is the critical point 
where the character of the order changes from first to 
second order. The critical point is shown by the black 
circle in Fig. 6. 

The change from second- to first-order transitions can 
be qualitatively explained by the Landau theory. 22 ) The 



free energy is then given by 



b i + c6 6 , 



(33) 



where 4> is an order parameter, and the coefficient c must 
be positive. The coefficient r is given by r = a(T — T c ) 
where T c is a second-order transition temperature. The 
character of the transition is of second order for b > 0, 
but becomes of first order for b < 0. The case with 6 = 
corresponds to the critical point. The Landau theory also 
explains the temperature dependence of order parame- 
ters shown in Fig. 7. At the critical point with b = 0, 
(j) near the transition point behaves as <f> ~ r 1 / 4 , where 
the critical exponent is 1/4 instead of 1/2. Therefore if 
b approaches as b — > +0, the initial rise becomes sharper 
reflecting the change of the exponent as in Fig. 7. 

5. Correlation Functions and Density of States 

5.1 Local Correlation Functions 

In this section, we elucidate the nature of the ordered 
state from correlation functions. Using pseudo spins, the 
local magnetic susceptibilities are defined by 



7(5 

Xm = 



[<T r S*(T)Sf> - {S${Si)] dr. (34) 



Since the magnetic dipole is given by J z — ^ 7 a 7 S* with 
appropriate a 7 , the magnetic susceptibility is represented 
as x.J — a -y a sXM- The parameter a 1 depends on the 
wave functions of the singlet-triplet states. 23 - 1 Let us dis- 
cuss the local magnetic susceptibility with A = shown 
in Fig. 11, which corresponds to a pair of the Kondo lat- 
tice models. In this case, the relations Xm = Xm anc ^ 
Xm = arc satisfied. Below the transition temperature 
T ~ 0.015, Xm splits into two values because of the emer- 
gence of the staggered CDW order. The B-sublattice has 
conduction electrons fewer than the A-sublattice. Then 
the B-sublattice shows a Curie-like behavior reflecting 
the localized spins with A = 0. The remaining entropy 
of the localized spins may lead to the magnetic insta- 
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Fig. 12. (color online) Local equal-time correlations (a) (Si • s c i) 
and (b) (Si ■ S 2 > with A = 0.2. 




Two Kondo 
Singlets CEF singlet 



Fig. 13. Schematic picture for the staggered order with Kondo 
and CEF singlets. The lower part represents localized electrons, 
while the upper part on the lines represents conduction electrons. 
The black and white circles show the presence and absence of 
electrons, respectively. 



bility at low temperatures. On the other hand, the sus- 
ceptibility for A-sublattice is strongly suppressed. This 
clearly indicates the formation of the Kondo singlet at 
A-sublattice. 

Next we consider the situation with the finite CEF 
splitting also shown in Fig. 11. In the case of A = 0.1, Xm 
is finite owing to the correlation between pseudo spins as 
shown in Fig. 11. The spin fluctuation at B-sublattice 
is suppressed, and all susceptibilities show paramagnetic 
behavior. 



Let us examine equal-time spin correlations that clar- 
ify properties of each sublattice with the finite CEF 
splitting. We first consider (Si ■ s c i), which is equal to 
(S2 -s C 2) under the present condition. Figure 12(a) shows 
the temperature dependence of (Si ■ s c i), which is en- 
hanced on A-sublattice and suppressed on B-sublattice. 
Bearing n c A > n C B in mind, we conclude that the local- 
ized spin on A-sublattice forms the Kondo singlet. On 
the other hand, the correlation (Si ■ S2) between local- 
ized spins is enhanced at B site as shown in Fig. 12(b). 
Hence, the B-sublattice corresponds to the CEF singlet. 

Let us estimate the magnitude of the effective CEF 
splitting for the CEF-singlet site. Taking the effective 
Hamiltonian for the localized states as T-L^ = KSi ■ S2, 
we obtain the susceptibility as Xm = — V(2A) for the 
ground state. Here A is the effective CEF splitting, which 
can be estimated in B-sublattice using the value in Fig. 
11. The result is A = 0.098, which is very close to the 
original CEF splitting A = 0.1. Hence, spatially ex- 
tended Kondo singlets do not significantly affect the mag- 
nitude of the CEF splitting. Besides, in Figure 12(b), 
(Si ■ S2) in the low-temperature limit is not far from 
—0.75 expected for the isolated singlet. Hence the CEF 
singlet is almost decoupled from conduction electrons. 
We note that the correspondence between A and A does 
not hold for A < 0.05. This corresponds to the fact that 
the CEF-singlet site tends to be magnetically polarized 
near A = 0. 

Thus, the present order with one conduction electron 
per site turns out to be a staggered order with the Kondo 
and CEF singlets. Figure 13 schematically shows this 
staggered order. Except for the strong coupling limit, 
the number of conduction electrons at the CEF singlet 
site is not zero because the Kondo singlets are spatially 
extended. 

5.2 Density of States 

The single-particle dynamics can be derived from the 
Green function (3) and (4). The density of states of con- 
duction electrons is given by 



pM 



A=A,B 



1 



ImGt(w + M) 



(35) 



We note that the Green functions does not depend on 
the labels 7 and a in the present condition. The Pade 
approximation is used for analytic continuation from 
imaginary Matsubara frequencies. 

Figure 14 shows the density of states for A = 0.2. In 
the disordered state (T = 0.04), the density of states 
shows the metallic behavior. The pseudo-gap structure 
around u> ~ 0.3 is interpreted as a kind of hybridization 
between conduction electrons and the localized states by 
the Kondo effect. Although there is no real hybridiza- 
tion because the pseudo spins do not have charge de- 
grees of freedom, strong renormalization by the Kondo 
effect gives rise to an electronic state that allows this 
interpretation. In the ordered state (T = 0.015), an en- 
ergy gap opens at the Fermi level. The sharp peak below 
the gap comes from the sublattice for the Kondo singlet, 
while the peak above the gap is due to the CEF singlet 
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Fig. 14. (color online) Density of states of conduction electrons 
in normal (T = 0.015) and ordered (T = 0.010) phases with 
A = 0.2. The bare density of states (J = 0) at quarter filling is 
also shown for comparison. 



site. Hence, this double-peak structure clearly shows the 
difference of the occupation number between the Kondo 
and CEF singlet sites. 

The origin of the insulating behavior is explained as 
follows. Each band has one conduction electron per unit 
cell in the ordered state. Provided that the localized spin 
at the Kondo-singlet site participates to the conduction 
band, each band is filled by two "electrons". Then the 
system can be an insulator as in the Kondo lattice at 
half filling. In this viewpoint, the staggered Kondo-CEF 
singlet order may be regarded as alternating itinerant 
and localized sites of /-electrons. 

6. Discussion 

6.1 Relation to CDW in Ordinary Kondo Lattice 

We discuss how the new electronic order found in the 
present paper is related to known orders in the ordi- 
nary Kondo lattice. In the limit of A = 0, the present 
2BSTKLM is reduced to a pair of Kondo lattice models. 
Let us briefly summarize the electronic order found in 
the ordinary Kondo lattice. In addition to the magnetic 
order, it has been found in ref. 13 that the Kondo lat- 
tice has a CDW order at quarter filling. In the strong 
coupling, the CDW is visualized as alternating Kondo 
singlets and localized spins. The corresponding electron 
number n c of conduction electrons per site is 1/2. As- 
sociated with hopping of conduction electrons, which is 
regarded as perturbation from the strong-coupling limit, 
a Kondo singlet and a local spin can exchange their po- 
sitions. This process is of first-order with respect to hop- 
ping. On the other hand, the second-order perturbation 
leads to inter-site attraction between a Kondo singlet 
and a local spin. 24,25 ' Although smaller attraction arises 
between Kondo singlets as well, the total second-order 
perturbation gives effective inter-site repulsion between 
Kondo singlets. Hence, this repulsive interaction gives a 
chance to stabilize the CDW order by partially sacrific- 
ing the hopping energy with intermediate coupling. In 
infinite dimensions, the CDW order is indeed stabilized 
at quarter filling with n c — 0.5 according to ref. 13. 



Since the uncompensated spin sites still have substan- 
tial entropy, the CDW alone cannot be the ground state 
in the Kondo lattice. It is likely that remaining spins form 
a magnetic order on top of the CDW background. 13 ) The 
magnetic fluctuation is clearly seen in the local magnetic 
susceptibility as shown in Fig. 11 with A = 0. Unless 
the CEF splitting is large enough, this situation may re- 
main in the 2BSTKLM with finite A. The region where 
we expect the magnetic ground state is roughly drawn in 
Fig. 6. However, we have not been able to demonstrate 
its existence, since the solution with small A and T in 
the two-sublattice DMFT does not converge. This indi- 
cates that the magnetic order has a longer periodicity 
than described by the two-sublattice system. 

Let us now consider the ordered phase in the 
2BSTKLM with larger A. We assume that the strong 
coupling limit in the 2BSTKLM is described by Kondo- 
singlet site and CEF-singlet site. In this case, the entropy 
vanishes even without magnetic order. One may regard 
the Kondo-singlet site as occupied by a fictitious spin- 
less fermion, and the CEF singlet site as vacant site of 
the fermion, i.e., a hole. In a similar manner to the or- 
dinary Kondo lattice in the strong coupling limit, the 
spinless fermions have an effective hopping and inter- 
site repulsion. This inter-site interaction tends to form a 
non-magnetic order, namely the staggered Kondo-CEF 
singlet order. It is clear that this non-magnetic order is 
a characteristic of non-Kramers systems. As shown in 
Figs. 6 and 8, moderate values of A stabilize the stag- 
gered Kondo-CEF singlet order. 

6.2 Itinerant and Localized Characters 

Next we discuss the staggered Kondo-CEF singlet or- 
der from the aspect of itinerant and localized characters 
of /-electrons. For f 1 system, the itinerant character is 
realized by the Kondo effect as heavy fermion state. If 
/ electrons are localized, on the contrary, a magnetic 
order appears by the RKKY interaction. The competi- 
tion between the Kondo effect and the RKKY interac- 
tion leads to the quantum phase transition between the 
ordered and disordered phases. For / 2 system with CEF 
singlet, both the itinerant and localized limits are disor- 
dered phases where the Kondo effect is dominant in the 
itinerant regime. The staggered Kondo-CEF singlet or- 
der is realized in the competing region, and interpreted 
as alternating sites of itinerant and localized states of /- 
electrons. In the weak-coupling, the itinerant character 
is responsible for the insulating ground state at quarter 
filling. 

The inter-site interaction leading to the staggered 
Kondo-CEF singlet order is different from the RKKY 
interaction. The effective repulsion between the Kondo 
singlets is the dominant mechanism for the present or- 
der. It is notable that the RKKY interaction is under- 
stood from the weak coupling limit, while the present 
staggered Kondo-CEF singlet order is understood natu- 
rally from the strong coupling limit. 

6.3 Relevance to Real Systems 

Let us finally discuss possible application of the 
present results to understanding PrFe4Pi2- This mate- 
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rial shows the Kondo-like behavior in the resistivity and 
undergoes a non-magnetic order at T = 6.5K. 26 ) In 
the ordered phase, a field-induced staggered moment is 
observed. 27 ) From phcnomenological and experimental 
analysis, this order is identified as a scalar order, 28 ~ 30 ) 
but the corresponding microscopic state is not yet clear. 
Inelastic neutron scattering experiment shows character- 
istic behaviors such as broad quasi-elastic peak in the 
disordered phase, and the inelastic peak in the ordered 
phase. 31 ~ 33 ) 

We remark that in the present staggered Kondo-CEF 
singlet order, the difference of the local susceptibility 
Xm between two sublattices results in an appearance 
of field- induced antiferromagnetic moment. Furthermore, 
our model naturally explains the appearance of CEF ex- 
citations only below the transition temperature. How- 
ever, the realistic band structure 34 ) is rather different 
from the identical two conduction bands taken in the 
present paper, which leads to an insulating ground state. 
More refinement is necessary for serious comparison with 
real systems, which will be given in separate publications. 

7. Summary and Outlook 

We have applied the DMFT combined with CT-QMC 
to the 2BSTKLM where CEF singlet-triplet states inter- 
act with two-band conduction electrons. The instability 
of the staggered ordered phase is derived by using the for- 
mulation of the susceptibility in two-sublattice systems. 
In the framework of the DMFT, physical quantities such 
as susceptibility, order parameter, correlation functions 
and the density of states have been calculated at finite 
temperatures. 

In the 2BSTKLM with one conduction electron per 
site, we have found the staggered order with Kondo and 
CEF singlets. The equal-time correlation shown in Fig. 
12 clearly shows this staggered ordering. This electronic 
order accompanies the CDW of conduction electrons be- 
cause they gather at the Kondo singlet site to screen 
the localized moments. Below the transition tempera- 
ture, the system becomes insulating as in the Kondo 
insulator, which is seen in the density of states. With 
different character of conduction bands, however, the in- 
sulating behavior should no longer hold. 

Although we have considered only magnetic and 
charge susceptibilities in this paper, it is also possible to 
calculate the pairing susceptibility for s-wave supercon- 
ductivity in the DMFT. In the 2BSTKLM, a pairing is 
possible mediated by the CEF excitation from singlet to 
triplet states. Indeed, we have observed in preliminary 
calculations an instability toward superconductivity at 
low temperatures. Here the singlet pairing between elec- 
trons in different conduction bands is realized. We shall 
discuss aspects related to the superconductivity in a sep- 
arate paper. 
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Appendix: Useful Formulae for Susceptibilities 

As we have seen in §2, the susceptibility in the lattice 
system is derived from the local susceptibility. In addi- 
tion to the calculation of the local susceptibility, we need 
to evaluate the susceptibilities without the vertex part. 
In this Appendix, we derive relevant formulae to calcu- 
late the local susceptibility. First of all, we define the 
following two complex functions: 

pO) g{z) 



Fi(z)= J . 



de- 



P(e) 



dg(z) 
dz : 



where 



«.) - / 



de 



(A-l) 
(A-2) 

(A-3) 



We have used the symmetric condition p(e) = p(—e). As 
shown later, the functions i*\ and F 2 are related to the 
staggered and uniform components, respectively. We can 
calculate the local Green function in the original Bril- 
louin zone from g(z). In the case with the hypercubic 
lattice, for example, g(z) is represented by an error func- 
tion. 14 ) 

The local Green function in the two-sublattice system 
is given by eq. (4). Only the diagonal elements survive 
the summation with respect to k, and are written as 



Gt(z) = (- x (z)F 1 (z), 



(A-4) 



with z = \/Ca(z)Cb(z). We omit the suffix a throughout 
this Appendix. On the other hand, the uniform suscep- 
tibility given in eq. (10) is calculated from the following 
form: 



Jjj2^2 G k A ( z i) G k A ( z 2) 
_(b(zi)(b{z 2 ) i 



z? z 2 



■[Fi(z 2 ) — Fi(%i) 



(A-5) 



— ^G k ( Zl )G k (z 2 ) = 

(A-6) 

In the special case with z\ = z 2 = z, we obtain 

(A-7) 

^E' GAB W G fe A W = -7^) - F ^)\ (A-8) 

These relations can also be obtained via z-derivative of 
eqs. (A-5) and (A-6). 

It is easy to confirm that these expressions reproduce 
the susceptibility in the original Brillouin zone. In the 
normal state, the relation (a = Cb = C is satisfied. Using 
eqs. (A-7) and (A-8), the uniform and staggered suscep- 
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tibility defined in eq. (6) is given by 

x u„if (z) = _p 2 (c(z)) ) (A . Q) 

X ^(z) = -F 1 (((z)). (A-10) 

Thus, the susceptibilities without vertex functions can 
be calculated from the functions i*\ and F 2 . 
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